A posteriori error covariances in variational data assimilation

The problem of variational data assimilation for a nonlinear evolution model is formulated as an optimal control problem to find some unknown parameters of the model. The equation for the error of the optimal solution is derived through the statistical errors of the input data (background, observation, and model errors). A numerical algorithm is developed to construct an a posteriori covariance operator of the analysis error using the Hessian of an auxiliary optimal control problem based on the tangent linear model constraints. The methods of data assimilation (DA) have become an important tool for analysis of complex physical phenomena in various fields of science and technology. These methods make it possible to combine mathematical models, data resulted from instrumental observations, and a priori information. The problems of variational data assimilation can be formulated as optimal control problems (e.g. [7, 9]) to find unknown model parameters such as initial and/or boundary conditions, righthand sides (forcing), and distributed coefficients. A necessary optimality condition reduces the problem to an optimality system which includes input errors; hence the error in the optimal solution. The statistical properties of the optimal solution error are important for estimating the efficiency of data assimilation in terms of reducing uncertainties in the model parameters and, therefore, in the model output. The error in the optimal solution can be derived through the errors in the input data using the Hessian of the cost functional of an auxiliary DA problem. For a deterministic case it has been done in [8]. If the errors in the input data are random and subjected to a normal distribution, then for a linearized problem (tangent linear approximation of the model) the covariance matrix of the analysis (optimal estimation of the initial condition) error is given by the inverse of the Hessian matrix of the cost functional (see e.g. [4, 5, 12, 14 –16]). This result was given (see e.g. [12]) for a discretized problem. In [3], a similar result was obtained for the continuous operator formulation. We showed that in the nonlinear case a posteriori covariance can often be approximated by the inverse Hessian of the auxiliary control problem (‘H-covariance’) beyond the validity of the tangent linear hypothesis (TLH). Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow 119333, Russia †MOISE project (CNRS, INRIA, UJF, INPG); LJK, Université Joseph Fourier, BP 51, 38051, Grenoble Cedex 9, France ‡Department of Civil Engineering, University of Strathclyde, John Anderson Building, 107 Rottenrow, Glasgow, G4 ONG, UK This paper presents a generalization of the theoretical results reported in [3] to the case of model errors. The equation for the error of the optimal solution is derived through the statistical errors of the input data (background, observation, and model errors). A numerical algorithm is developed to construct an a posteriori covariance operator of the analysis error using the Hessian of an auxiliary optimal control problem based on the tangent linear model constraints. Different approaches to model error formulation in 4D-Var are presented in [1, 17] (see also citations in [1]). The paper is organized as follows. In Section 1, we give the statement of the variational data assimilation problem for a nonlinear evolution model to identify the model parameters. In Section 2, the equation of the error of the optimal solution is derived through the errors of the input data. In Section 3, we derive the formulas and the algorithm for constructing the covariance operator of the optimal solution errors through the covariance operators of the input errors using the Hessian of the cost functional of the auxiliary control problem. 1. Statement of the problem Consider the mathematical model of a physical process that is described by the evolution problem ( ∂φ ∂ t = F(φ ;λ )+ f ; t 2 (0;T ) φ t=0 = u (1.1) where φ = φ(t) is the unknown function belonging for any t to a Hilbert space X , u 2 X , F is a nonlinear operator mapping Y Yp into Y with Y = L2(0;T ;X), k kY = ( ; )1=2 Y , Yp is a Hilbert space (the space of parameters, or control space), f 2Y . Suppose that for given u 2 X ; f 2Y and λ 2Yp there exists a unique solution φ 2Y to (1.1). The function λ is assumed to be an unknown parameter of the model. Let us introduce the functional S(λ ) = 1 2 (V1(λ λb);λ λb)Yp + 12 (V2(Cφ φobs);Cφ φobs)Yobs (1.2) where λb 2Yp is a prior (background) function, φobs 2Yobs is a prescribed function (observational data), Yobs is a Hilbert space (observation space), C : Y ! Yobs is a linear bounded operator, V1 : Yp ! Yp and V2 : Yobs ! Yobs are symmetric positive definite operators. Consider the following data assimilation problem with the aim to identify the parameter λ : find λ 2 Yp and φ 2 Y such that they satisfy (1.1), and the functional A posteriori error covariances 163 S(λ ) on the set of solutions to (1.1) takes the minimum value, i.e. 8>>><>>>: ∂φ ∂ t = F(φ ;λ )+ f ; t 2 (0;T ) φ t=0 = u S(λ ) = inf v2Yp S(v): (1.3) We suppose that the solution of (1.3) exists. (Note that the solvability of optimal control parameter estimation problems has been addressed, e.g., in [2, 11].) The necessary optimality condition S0(λ ) = 0 reduces problem (1.3) to the optimality system [13]: 8<: ∂φ ∂ t = F(φ ;λ )+ f ; t 2 (0;T ) φ t=0 = u (1.4) ( ∂φ ∂ t (F 0 φ(φ ;λ )) φ = C V2(Cφ φobs); t 2 (0;T ) φ t=T = 0 (1.5) V1(λ λb) (F 0 λ (φ ;λ )) φ = 0: (1.6) Here F 0 φ(φ ;λ ) : Y !Y and F 0 λ (φ ;λ ) : Yp!Y are the Frechet derivatives of F with respect to φ and λ , respectively, (F 0 φ(φ ;λ )) :Y !Y; (F 0 λ (φ ;λ )) :Y !Yp are their adjoints, and C is the adjoint to C defined by (Cφ ;ψ)Yobs = (φ ;C ψ)Y ; φ 2Y;ψ 2 Yobs. We assume that system (1.4)–(1.6) has a unique solution. Suppose that λb = λ̄ +ξ1; φobs =Cφ̄ +ξ2, f = f̄ +ξ3, where ξ1 2 Yp; ξ2 2 Yobs, ξ3 2 Y , and φ̄ is the (‘true’) solution to the problem (1.1) with λ = λ̄ and f = f̄ : 8<: ∂ φ̄ ∂ t = F(φ̄; λ̄ )+ f̄ ; t 2 (0;T ) φ̄ t=0 = ū: (1.7) The functions ξ1;ξ2;ξ3 are treated as the errors of the input data λb;φobs; f (‘background’, observation, and model errors, respectively). For V1 and V2 in (1.2), one usually has V1 = V 1 ξ1 ; V2 = V 1 ξ2 , where Vξ is the covariance operator of the corresponding error ξ . 2. Equation for the optimal solution error Let us derive the equation for the optimal solution error through input errors. Let δφ = φ φ̄; δλ = λ λ̄ . Then, from (1.7) and the optimality system (1.4)–(1.6), we obtain ∂δφ ∂ t Fφ φ̃ λ̃ δφ Fλ φ̃ λ̃ δλ ξ3 t 0;T ) δφ t 0 0 (2.1) ∂φ ∂ t Fφ φ λ φ C V2 Cδφ ξ2 t 0;T ) φ t=T = 0 (2.2) V1(δλ ξ1) (F 0 λ (φ ;λ )) φ = 0 (2.3) where φ̃ = φ̄ + τ(φ φ̄); λ̃ = λ̄ + τ(λ λ̄)τ 2 [0;1℄: Note that φ̃ = φ̄ + τδφ ; φ = φ̄ +δφ , λ̃ = λ̄ + τδλ ; λ = λ̄ +δλ . The system (2.1)–(2.3) can be written in the form: ( ∂δφ ∂ t F 0 φ(φ̄ ; λ̄ )δφ = F 0 λ (φ̄ ; λ̄ )δλ +ξ3 + ξ̃3; t 2 (0;T ) δφ jt=0 = 0 (2.4) ( ∂φ ∂ t (F 0 φ(φ̄ ; λ̄ )) φ = C V2(Cδφ ξ2)+ξ4; t 2 (0;T ) φ t=T = 0 (2.5) V1(δλ ξ1) (F 0 λ (φ̄ ; λ̄ )) φ = ξ5 (2.6) where ξ̃3 = [F 0 φ(φ̃; λ̃ ) F 0 φ(φ̄ ; λ̄ )℄δφ +[F 0 λ (φ̃ ; λ̃ ) F 0 λ (φ̄; λ̄ )℄δλ ξ4 = [(F 0 φ(φ ;λ )) (F 0 φ(φ̄ ; λ̄ )) ℄φ ; ξ5 = [(F 0 λ (φ ;λ )) (F 0 λ (φ̄ ; λ̄ )) ℄φ : Let us introduce the operator H :Yp!Yp defined by the successive solutions of the following problems: ( ∂ψ ∂ t F 0 φ(φ̄ ; λ̄ )ψ = F 0 λ (φ̄ ; λ̄ )v; t 2 (0;T ) ψ jt=0 = 0 (2.7) ( ∂ψ ∂ t (F 0 φ(φ̄ ; λ̄ )) ψ = C V2Cψ ; t 2 (0;T ) ψ t=T = 0 (2.8) Hv=V1v (F 0 λ (φ̄; λ̄ )) ψ : (2.9) Below we introduce four auxiliary operators R1;R2;R3;R4. Let R1 = V1. Let us introduce the operator R2 :Yobs !Yp acting on the functions g 2Yobs according to the formula R2g= (F 0 λ (φ̄ ; λ̄ )) θ (2.10) A posteriori error covariances 165 where θ is the solution to the adjoint problem ( ∂θ ∂ t (F 0 φ(φ̄ ; λ̄ )) θ = C V2g; t 2 (0;T ) θ t=T = 0: (2.11) The operator R3 :Y !Yp is defined on the functions q 2Y as follows: ( ∂θ1 ∂ t F 0 φ(φ̄ ; λ̄ )θ1 = q; t 2 (0;T ) θ1jt=0 = 0 (2.12) ( ∂θ 1 ∂ t (F 0 φ(φ̄ ; λ̄ )) θ 1 = C V2Cθ1; t 2 (0;T ) θ 1 t=T = 0 (2.13) R3q= (F 0 λ (φ̄; λ̄ )) θ 1 : (2.14) The operator R4 :Y !Yp is defined on the functions h 2Y as ( ∂θ 2 ∂ t (F 0 φ(φ̄ ; λ̄ )) θ 2 = h; t 2 (0;T ) θ 1 t=T = 0 (2.15) R4h= (F 0 λ (φ̄; λ̄ )) θ 2 : (2.16) From (2.7)–(2.16) we conclude that system (2.4)–(2.6) is equivalent to the single equation for δλ : Hδλ = R1ξ1+R2ξ2+R3(ξ3+ ξ̃3)+R4ξ4+ξ5: (2.17) This is the exact equation for δλ . Under the hypothesis that H is invertible, we get δλ = T1ξ1+T2ξ2+T3(ξ3+ ξ̃3)+T4ξ4+T5ξ5 (2.18) where Ti=H 1Ri; i= 1;2;3;4; T5 =H 1; T1 :Yp!Yp; T2 :Yobs!Yp; T3;T4 :Y!Yp: Note, however, that the functions φ ;λ ; φ̃; λ̃ in (2.1)–(2.3) depend on ξ1;ξ2;ξ3, as the result, the terms T3ξ3;T4ξ4;T5ξ5 also depend on ξ1;ξ2;ξ3 (nonlinearly), and it is not possible to represent δλ through ξ1;ξ2;ξ3 in an explicit form. To derive from (2.18) the covariance operator of δλ , we need to introduce some approximation of (2.18). Since φ̃ = φ̄ + τδφ ; φ = φ̄ + δφ , λ̃ = λ̄ + τδλ ; λ = λ̄ + δλ , we assume that T3ξ̃3 0; T4ξ4 0; T5ξ5 0 (2.19) then (2.18) reduces to δλ = T1ξ1+T2ξ2+T3ξ3 (2.20) which is equivalent to the system: ( ∂δφ ∂ t F 0 φ(φ̄ ; λ̄ )δφ = F 0 λ (φ̄ ; λ̄ )δλ +ξ3; t 2 (0;T ) δφ jt=0 = 0 (2.21) ( ∂φ ∂ t (F 0 φ(φ̄ ; λ̄ )) φ = C V2(Cδφ ξ2); t 2 (0;T ) φ t=T = 0 (2.22) V1(δλ ξ1) (F 0 λ (φ̄ ; λ̄ )) φ = 0: (2.23) Taking into account the definition of ξ̃3;ξ4;ξ5; it is easily seen that assumption (2.19) comes from the first-order approximation of the Taylor–Lagrange formula under the hypothesis that F is twice continuously Frechet differentiable [10]. Using this formula, the errors ξ̃3;ξ4;ξ5; may be expressed through the second derivatives of F , and the values of the norms of T3ξ̃3;T4ξ4;T5ξ5 can be estimated, thus giving the possibility to analyse the approximation error when taking (2.20) instead of (2.18). The problem (2.21)–(2.23) is a linear data assimilation problem; for fixed λ̄ ; φ̄ it is the necessary optimality condition to the following (auxiliary) minimization problem: Find δλ and δφ such that 8>><>>: ∂δφ ∂ t F 0 φ(φ̄ ; λ̄ )δφ = F 0 λ (φ̄ ; λ̄ )δλ +ξ3; t 2 (0;T ) δφ jt=0 = 0 S1(δλ ) = inf v2YpS1(v) (2.24) where S1(δλ ) = 1 2 (V1(δλ ξ1);δλ ξ1)Yp + 12 (V2(Cδφ ξ2);Cδφ ξ2)Yobs : (2.25) The Hessian H of functional (2.25) is defined on v2Yp by (2.7)–(2.9). Note that for ξ2 = 0 the operator H coincides with the HessianH of the original nonlinear DA problem on the exact solution λ̄ . The Hessian H acts in Yp as a self-adjoint operator with the domain of definition D(H) =Yp. Moreover, due toV1;V2, the operator H is positive definite, and hence invertible. Note that if the tangent linear hypothesis is valid (e.g. [4]), then for small δφ ;δλ we can choose (2.19). However, the transition from (2.18) to (2.20) may not necessarily require the tangent linear hypothesis to be valid. As follows from (2.20), the influence of the errors ξ1;ξ2;ξ3 on the value of the error δλ of the optimal solution is determined by the operators H 1R1;H 1R2; H R3, respectively. The values of the norms of these operators can be considered as sensitivity coefficients: the less is the norm of the operator H Ri, the less impact on δλ is given by the corresponding error ξi. This criterion was used for deterministic error analysis in [6, 8] with the aim to identify the initial condition. Here, assuming the statistical structure of the errors ξ1;ξ2;ξ3, we will derive the covariance operator of the optimal solution (parameter) error through the covariance operators of the input errors and develop a numerical algorithm to construct the covariance operator of the optimal solution error using the covariance operators of the input errors. 3. Covariance operator as the inverse Hessian Consider error equation (2.20), where Ti = H 1Ri; i = 1;2;3; T1 : Yp!Yp; T2 : Yobs!Yp; T3 : Y!Yp: Below we suppose that the errors ξ1;ξ2;ξ3 are normally distributed, unbiased, and mutually uncorrelated. By Vξi we denote the covariance operator of the corresponding error ξi; i= 1;2;3, i.e. Vξ1 = E[( ;ξ1)Ypξ1℄; Vξ2 = E[( ;ξ2)Yobsξ2℄; Vξ3 = E[( ;ξ3)Yξ3℄, where E is the expectation. By Vδλ we denote the covariance operator of the optimal solution (analysis) error: Vδλ = E[( ;δλ )Ypδλ ℄. From (2.20) we get Vδλ = T1Vξ1T 1 +T2Vξ2T 2 +T3Vξ3T 3 : (3.1) To find the covariance operator Vδλ , we need to construct the operators TiVξiT i ; i= 1;2;3. We proved in [13] that T1Vξ1T 1 +T2Vξ2T 2 = H 1 (3.2) where H is the Hessian of the auxiliary data assimilation problem (2.24)–(2.25) defined by (2.7)–(2.9). Then, Vδλ = H 1 +T3Vξ3T 3 : (3.3) Consider now the operator T3 =H R3. To construct T3Vξ3T 3 , we need to derive R 3. For q 2Y; p 2Yp, we have from (2.12)–(2.14): (R3q; p)Yp = ((F 0 λ (φ̄ ; λ̄ )) θ 1 ; p)Yp = (C V2Cθ1;φ)Y = (V2Cθ1;Cφ)Yobs where θ1;θ 1 are the solutions to (2.12)–(2.13), and φ is the solution to (2.7) for v= p. Further, (R3q; p)Yp = (θ1;C V2Cφ)Y = (q;φ )Y and R 3p= φ , where φ is the solution to the adjoint problem: ( ∂φ ∂ t (F 0 φ(φ̄ ; λ̄ )) φ = C V2Cφ ; t 2 (0;T ) φ t=T = 0: (3.4) Let Q = R3Vξ3R 3 : Yp ! Yp. The operator Q can be defined as follows: for a given p 2 Yp find φ as the solution of (2.7) for v = p, find φ as the solution of (3.4), and for q =Vξ3 φ find θ1;θ 1 as the solutions of (2.12)–(2.13); then put R3Vξ3R 3 = (F 0 λ (φ̄; λ̄ )) θ 1 . Therefore, T3Vξ3T 3 = H 1R3Vξ3R 3H 1 = H QH , and Q is defined by the successive solutions of the following problems (for the given p 2Yp): ( ∂φ ∂ t F 0 φ(φ̄; λ̄ )φ = F 0 λ (φ̄; λ̄ )p; t 2 (0;T ) φ jt=0 = 0 (3.5) ( ∂φ ∂ t (F 0 φ(φ̄ ; λ̄ )) φ = C V2Cφ ; t 2 (0;T ) φ t=T = 0 (3.6) ( ∂θ1 ∂ t F 0 φ(φ̄; λ̄ )θ1 = Vξ3 φ ; t 2 (0;T ) θ1jt=0 = 0 (3.7) ( ∂θ 1 ∂ t (F 0 φ(φ̄ ; λ̄ )) θ 1 = C V2Cθ1; t 2 (0;T ) θ 1 t=T = 0 (3.8) then Qp= (F 0 λ (φ̄ ; λ̄ )) θ 1 : (3.9) The algorithm (3.5)–(3.9) can be used to compute the operator Q numerically. Then, from (3.3), we come to the main result of the paper. Theorem 3.1. The covariance operator Vδλ of the optimal solution error is given by the formula Vδλ = H 1+H QH 1 (3.10) where H is the Hessian of the functional S1 defined by (2.7)–(2.9), and the operator Q is defined by (3.5)–(3.9). It is not difficult to show that Q is the Hessian of the following minimization problem: Find v 2Yp such that SQ(v) = inf p2YpSQ(p) (3.11) where SQ(p) = 1 2 (Vξ3 φ ;φ )Y (3.12) and φ is defined through p by the successive solutions of the following problems: ( ∂φ ∂ t F 0 φ(φ̄; λ̄ )φ = F 0 λ (φ̄; λ̄ )p; t 2 (0;T ) φ jt=0 = 0 (3.13) ( ∂φ ∂ t (F 0 φ(φ̄ ; λ̄ )) φ = C V2Cφ ; t 2 (0;T ) φ t=T = 0: (3.14)